namespace hms::swe {

template<typename MeshType>
template<BlockLocation loc>
void SolverThread<MeshType>::addMassSources(
	const Block<loc>& block,
	Array3XsMap& qMap
){
	/* mass sources can be added to q directly, because the preconditioning for
	 * friction only applies to the momentum terms.
	 * However, for implementing momentum sources, the sources must be added
	 * to the flux vector instead. */
	scalar dt { solver().time().step() };
	bool hasMassSource {false};

	solver().visitUsedSourceFields( [&]( const auto& sourceField ){
		hasMassSource = true;
		block.forEachRow( [&](Index y){
			auto qRow { block.rowMap_mem(qMap, y).template topRows<1>() };
			auto sourceRow = [&](){
				if constexpr ( remove_qualifiers<decltype(sourceField)>::isUniform ){
					return
						block.template writeBlock<1>( sourceField.values() )
						.col(y).transpose();
				} else {
					return block.template rowMap<1>(
						block.template writeBlock<1>( sourceField.values() ), y
					);
				}
			}();
			qRow += dt * sourceRow;
		} );
	} );
	/* prevent negative water depths (may arise from infiltration) */
	if (hasMassSource){
		block.forEachRow( [&](Index y){
			auto qRow { block.rowMap_mem(qMap, y).template topRows<1>() };
			qRow = max(qRow, 0);
		} );
	}

	if constexpr ( !(
		std::is_same_v<MeshType, UniMesh> ||
		std::is_same_v<MeshType, RectMesh>
	) ){
		assert(false);
	}
}

} // namespace hms::swe